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Abstract 

Network reconstruction, i.e., obtaining network structure from data, is a central theme in systems bi- 
ology, economics and engineering. In some previous work, we introduced dynamical structure functions 
as a tool for posing and solving the problem of network reconstruction between measured states. While 
recovering the network structure between hidden states is not possible since they are not measured, in 
many situations it is important to estimate the minimal number of hidden states in order to understand 
the complexity of the network under investigation and help identify potential targets for measurements. 
Estimating the minimal number of hidden states is also crucial to obtain the simplest state-space model 
that captures the network structure and is coherent with the measured data. This paper characterizes 
minimal order state-space realizations that are consistent with a given dynamical structure function by 
exploring properties of dynamical structure functions and developing an algorithm to explicitly obtain 
such a minimal realization. 

I. Introduction 

Networks have received an increasing amount of attention in the last decade. In our "information- 
rich" world, the questions of network reconstruction and network analysis become crucial for 
the understanding of complex systems such as biological, social, or economical networks. In 
particular, the analysis of molecular networks has gained significant interest due to the recent 
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explosion of publicly available high-throughput biological data. In this context, the question of 
identifying and analyzing the network structure at the origin of measured data becomes a key 
issue. 

In some occasions, measured data is given in the form of input-output time-series that describes 
the effect of inputs on outputs (measured states) of a network. When data is generated by 
a linear system, a matrix transfer function describing the dynamic input-output behavior is 
generally obtained using system identification [1]. If the original state-space model is available 
or deducible, then the associated network structure can be readily obtained from it. However, 
a transfer function cannot, in general, recover, or realize, the original state-space model since 
the realization problem does not typically have a unique solution, i.e., different state-space 
realizations can generate the same input-output behavior. Since each of these realizations may 
suggest entirely different network structures, it is in general impossible to identify network 
structures from transfer functions alone. Therefore, more information, beyond input-output data 
used to identify a transfer function, is needed to prefer one state- space realization over another 
as a description of a particular system (2). 

Another difficulty in the network reconstruction problem comes from the fact that the real- 
ization problem becomes ill posed when some of the states are unobservable or "hidden" (this 
happens with just one hidden state |3} pp. 78]). As a result, failure to explicitly acknowledge 
the presence of hidden states and the resulting ambiguity in network structures can lead to 
a deceptive and erroneous process for network structure discovery. Consequently, determining 
from measured data the presence or absence of a causal relationship between two variables in a 
network is a challenging question. 

Motivated by this, we are focusing on the effect of hidden states in the network that we 
are aiming to reconstruct. A new representation for LTI systems, called dynamical structure 
functions was introduced in [4] and developed in (5j-|10) . Dynamical structure functions capture 
information at an intermediate level between transfer function and state space representation (see 
Figure [T]). Specifically, dynamical structure functions not only encode structural information at 
the measurement level, but also contain some information about hidden states. Based on the 
theoretical results presented in (4), we proposed some guidelines for the design of an experimental 
data- acquisition protocol which allows the collection of data containing sufficient information for 
the network structure reconstruction problem to become solvable. In particular, we have shown 
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that if nothing is known about the network, then reconstruction is impossible. If, however, 
one can make the following reasonable assumptions about the data-collection experiments: 
Networks have received an increasing amount of attention in the last decade. In our "information- 
rich" world, the questions of network reconstruction and network analysis become crucial for 
the understanding of complex systems such as biological, social, or economical networks. In 
particular, the analysis of molecular networks has gained significant interest due to the recent 
explosion of publicly available high-throughput biological data. In this context, the question of 
identifying and analyzing the network structure at the origin of measured data becomes a key 
issue. 

In some occasions, measured data is given in the form of input-output time-series that describes 
the effect of inputs on outputs (measured states) of a network. When data is generated by 
a linear system, a matrix transfer function describing the dynamic input-output behavior is 
generally obtained using system identification JTJ. If the original state-space model is available 
or deducible, then the associated network structure can be readily obtained from it. However, 
a transfer function cannot, in general, recover, or realize, the original state-space model since 
the realization problem does not typically have a unique solution, i.e., different state-space 
realizations can generate the same input-output behavior. Since each of these realizations may 
suggest entirely different network structures, it is in general impossible to identify network 
structures from transfer functions alone. Therefore, more information, beyond input-output data 
used to identify a transfer function, is needed to prefer one state- space realization over another 
as a description of a particular system Q. 

Another difficulty in the network reconstruction problem comes from the fact that the real- 
ization problem becomes ill posed when some of the states are unobservable or "hidden" (this 
happens with just one hidden state |3} pp. 78]). As a result, failure to explicitly acknowledge 
the presence of hidden states and the resulting ambiguity in network structures can lead to 
a deceptive and erroneous process for network structure discovery. Consequently, determining 
from measured data the presence or absence of a causal relationship between two variables in a 
network is a challenging question. 

Motivated by this, we are focusing on the effect of hidden states in the network that we 
are aiming to reconstruct. A new representation for LTI systems, called dynamical structure 
functions was introduced in [4] and developed in Q-|10) . Dynamical structure functions capture 
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Fig. 1. Mathematical structure of the network reconstruction problem using dynamical structure functions. Red arrows mean 
"uniquely determine", blue arrows indicate our work. 

information at an intermediate level between transfer function and state space representation (see 
Figure [T]). Specifically, dynamical structure functions not only encode structural information at 
the measurement level, but also contain some information about hidden states. Based on the 
theoretical results presented in Q, we proposed some guidelines for the design of an experimental 
data- acquisition protocol which allows the collection of data containing sufficient information for 
the network structure reconstruction problem to become solvable. In particular, we have shown 
that if nothing is known about the network, then reconstruction is impossible. If, however, one 
can make the following reasonable assumptions about the data-collection experiments: 
(A.l) for a network composed of p measured species, the same number of experiments p must 
be performed; 

(A.2) each experiment must independently control a measured species, i.e., control input i must 

first affect measured species z, 
then reconstruction is possible. Moreover, failure to meet the necessary informativity conditions 
results in a situation where any internal network structure fits the data equally well (e.g. a fully 
decoupled network or a fully connected network). If biologists have already some information 
about the network, as it is usually the case, then these conditions can be relaxed as explained 
in Q. Using dynamical structure functions as a mean to solve the network reconstruction 
problem, the following aspects need to be considered (see Figure [TJ: 

First (see (A) in Figure [T]), the properties of a dynamical structure function and its relationship 
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with the transfer function associated with the same system need to be precisely established (this 
was done in Q). 

Second (see (B) in Figure [Tf, an efficient method is developed to reconstruct networks in the 
presence of noise and nonlinearities (this was done in (TTJ). This method relies on the assumption 
that the conditions for network reconstruction presented above in (A.l) and (A.2) have been met. 
In our approach, we use the same information as traditional system identification methods, i.e., 
input-output data. However, with our method, steady-state (resp. time-series data) can be used 
to reconstruct the Boolean (resp. dynamical network) structure of the system (see jTTJ for more 
details). 

Third (see (C) in Figure [T]), once the dynamical structure function is obtained, as a main 
result of this paper, an algorithm for constructing a minimal order state- space representation 
consistent with such function is developed. In an application, this provides a way to estimate the 
complexity of the system by determining the minimal number of hidden states in the system. 
For example, in the context of biology it helps understand the number of unmeasured molecules 
in a particular pathway: a low number means that most molecules in that pathway have been 
identified and measured, showing a good understanding of the system; while a large number 
shows that there are still many unmeasured variables, suggesting that new experiments should 
be carried out to better characterize that pathway. 

The outline of the paper is as follows. Section |TT] reviews the definition of dynamical structure 
functions and their properties. The main result can be found in Section [III] where we propose 
a minimal order realization algorithm based on state-space realizations and pole-zero analysis. 
Simulation and discussion are addressed in Section [TV} Finally conclusions are presented in 
Section [V] 

II. System Model 

Consider a linear system (it can also be a linearization of some original nonlinear system) 
x = Ax + Bu, y = Cx. The transfer function associated with this system is given by G(s) = 
C(sl— A)~ l B. Typically, we can use standard system identification tools [ 1 ] to identify a transfer 
function G{s) from input-output data. 

Like system realization, network reconstruction also begins with the identification of a transfer 
function, but it additionally attempts to determine the network structure between measured states 
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without imposing any additional structure on the hidden states. As we have shown in [4], this 
requires a new representation of linear time-invariant systems: the dynamical structure function 
(defined later). An algorithm allowing the dynamical structure function to be obtained from 



input-output data is proposed in jTTJ. This paper assumes that the dynamical structure function 
has already been obtained from data, and will focus on finding one of its minimal state-space 
realizations. 

The dynamical structure function is obtained as follows: First, we transform [A,B,C] to 
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without changing G(s), where p = rank(C). The linear system dynamics then 
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where x = (y, z) G M n ° is the full state vector, y G MP is a partial measurement of the state, 
z are the n° — p "hidden" states, and u G R m is the control input. In this work we restrict 
our attention to situations where output measurements constitute partial state information, i.e., 
p < n°. We consider only systems with full rank transfer functions that do not have entire rows 
or columns of zeros, since such "disconnected" systems are somewhat pathological and only 
serve to complicate the exposition without fundamentally altering our conclusions. 
Taking the Laplace transforms of the signals in ([T]) yields 
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where Y, Z, and U are the Laplace transforms of y, z, and u, respectively. Solving for Z gives 

Z=(sl- A*,)- 1 A° 21 Y + (si - A° 22 y l B° 2 U 
Substituting this last expression of Z into ([2]) then yields 

sY = W°Y + V°U (3) 

where W° = A° n + A° 12 (si - A^)' 1 A° 21 and V° = B\ + A° 12 (si - A° 22 )~ x B° 2 . 

Now, let R° be a diagonal matrix formed of the diagonal terms of W° on its diagonal, i.e., 
R° = diag{V^ } = diag(WJi, W 22 , W° p ). Subtracting R°Y from both sides of ©, we obtain: 

(si -R°)Y = (W° -R°)Y + V°U 
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Note that W° — R° is a matrix with zeros on its diagonal. We thus have: 

Y = QY + PU (4) 

where 

Q = (si - R°Y l (W° - R°) (5) 

and 

P= (sI-R^'W (6) 

Note that Q is zero on the diagonal. 

Definition 1: Given the system ([T]), we define the dynamical structure function of the system 
to be [Q,P]. 

Note that, in general, Q(s) and P(s) carry a lot more information than G(s). This can be 
seen from the equality G(s) = (/ — Q(s))~ 1 P(s) (see [4] for details). However, Q(s) and P(s) 
carry less information than the state-space model ([T]) (see |TTJ). 

Definition 2: A dynamical structure function, [Q, P], is said to be consistent with a particular 
transfer function, G, if there exists a realization of G, of some order, and of the form ([T]), such 
that [Q, P] are specified by ([5]) and ([6]). Likewise, a realization is consistent with [Q, P] if that 
realization gives [Q, P] from ([5]) and ([6]). 

Definition 3: We say that a realization is G minimal if this realization corresponds to a 
minimal realization of G. We say that a realization is [Q,P] minimal if this realization is 
consistent with [Q, P] and its order is smaller than or equal to that of all realizations consistent 
with [Q,P]. 

The underlying principle to find a [Q, P] minimal realization is to search for a realization 
with the minimal number of hidden states. Such a realization is characterized by the minimal 
number of pole-zero cancellations in the transfer functions Q and P. 

Proposition 1: Given a dynamical system ([T]) and the associated dynamical structure functions 
[Q, P] with P° constructed as explained above (see ([T])-([6])), the following conditions must hold 

diag^J = lim (7) 

A° 11 -diag{A° u }= limsQ(s); (8) 

P x ° = lim sP(s). (9) 
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Proof: Eq. ([7]) is directly obtained from the definition of R°(s)\ 
lim R°{s) = lim dmg{W°(s)} 

s— >oo s— >oo 

= diag{ lim W°(s)} = diag{^J 

Since the proofs for eq. ([8]) and ([9]) are very similar, we focus on eq. ([8]) only. Using the fact 
that for any square matrix M, if M n — )► when n —t +00, then (/ — M) _1 YliLo M\ we 
obtain, from the definition of Q given in ([5]), Q(s) = J2Zi s ~ iR ° (W°(s) - R°(s)) and 

W°(s) = A° n + J2Zi s~ % A\ 2 A° 2 ^- Y A° 2Y , when s +oc. Hence, Q(s) = (A° u - R^s^s' 1 +r(s), 
in which r(s) is a matrix polynomial of 5, whose largest degree is —2. Finally, multiplying by 
s on both sides and taking the limit as s goes to 00 results in eq. ([8]). A similar argument can 
be used to prove eq. ([9]). ■ 

Remark 1: This Proposition reveals an important property of dynamical structure functions: 
they encode the direct causal relationships between observed variables. 

We present hereafter an illustrative example to help fix the ideas. 

Example 1: Consider a network with the structure depicted in Fig. [2j A linear state-space 
representation of this network is given by 
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where I 3 is the 3 x 3 identity matrix. Following the definitions in ([5]) and ([6]), we can write down 
the corresponding dynamical structure function [Q,P] as 
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Consistent with Proposition [TJ we can check using the above that: 
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Generally, there exist many realizations consistent with [Q, P]. In the following section, 



I 



^, i.e., a realization which is 



we focus on finding a [Q, P] minimal realization ^4, P, 
consistent with [Q, P] and which has minimal order, that is, with minimal dimension for A (and 
hence the lowest possible complexity). 



III. Algorithm to find a [Q, P] minimal realization 

From a dynamical structure function [Q, P] we cannot reconstruct [V^°,y°] since there is 
no information regarding the diagonal transfer function matrix R°. Given [Q, P] and a diagonal 
proper transfer function matrix R, a minimal realization of [W V] = [(sI — R)Q + R (sI — R)P] 
can be obtained as follows. We know that: 



[W V\ = [A n B 1 ]+A 12 (sI-A 22 )~ 1 [A 21 B 2 ] 



(10) 
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Fig. 2. (a) An example system with two inputs, three measured states (states 1, 2, and 3) and two hidden states (states 4 and 
5). (b) The corresponding dynamical network structure. 
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The idea is to start with an arbitrarily chosen P, and then use a state-space realization approach 
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Lemma 1: Suppose W,V,A= " and B = L are defined as in eq. ( [TO] ), then 

V and G share the same zeros. 

Proof: Since 5/ — W is the Schur complement of si — A 22 in si — A, then 

^{sI-W)^ ^-f\ . (11) 
v ; det(sI-A 22 ) 

Recall that V = B 1 + A 12 (sl - A 22 )~ 1 B 2 . Since (si - W)G = V, we thus have that V and G 

share the same zeros (15} page 153]. ■ 

Given a dynamical structure function [Q, P], a random choice of a proper diagonal transfer 

function matrix R is likely to result in additional zeros in V = (si - P)P. From Lemma [TJ 

this will lead to additional zeros in G which are associated to uncontrollable eigenvalues of the 

considered realization fl2\ Section 4] and of course does not lead to a minimal realization in 

eq. fLO] ). At this stage the following question arises: how can we find a proper diagonal transfer 

function matrix P* such that a minimal realization of [W V] is a [Q, P] minimal realization, 

i.e., 

P* =argmin^ deg {(s/ - P)^ 1 ^ sP] + [R 0]}, (12) 

where deg is the McMillan degree (3). Note that, since there are many choices for P* that 
minimize the order of minimal realizations of [W V], a chosen P* may be different from R°. 

Assume that all elements in [Q P] only have simple poles. This assumption can be relaxed 
but we adopt it here for simplicity. Also assume that [Q P] does not possess any poles at 
(otherwise we can change eq. ( fT2] ) to (si — R)(s — a) -1 [(s — a)Q (s — a)P] + [P 0], where 
a e R is not a pole of [Q P]). 

Proposition 2: Assume [I — Q P] only has simple poles and does not have any zero^} A 
minimal order realization of [W V] in ( [T0| ) can be achieved using a constant diagonal matrix 
P*. 

Proof: Assume P* has at least one term on the diagonal with the degree of the numerator 



greater or equal to 1, e.g., suppose the i th term in (si — P*)s 



-1 _ (s+b)ej(s) 
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These assumptions can be relaxed, see Section 3.6 of [14] for more details. 
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and deg{ei{s)) = deg((j>i(s)) > 1, where deg(-) returns the degree of a polynomial. Hence, the 
product (si — sP] will introduce deg(cj)i(s)) new poles and, due the assumption of 

simple poles, can at most eliminate deg(ei(s)) = deg((/>i(s)) poles since [I — Q P] does not 
have any zeros. As a consequence, we can change the i th term from ^^y^ to ^ without 
increasing the order. Doing this along all the elements of R* proves the result. ■ 
If R* is a constant matrix, the term [R* 0] in eq. ( fT2] ) is also a constant matrix. Therefore, 
the order of a minimal realization is only determined by (si — R*)s~ 1 [sQ sP] = N[sQ sP]. 
Thus, finding the "optimal" R* which leads to the minimal order in eq. ([12]) is equivalent to 
finding a diagonal proper transfer matrix TV [^] such that N[sQ sP] has as few poles as possible. 
Based on this idea, the following algorithm is proposed: 
Step 1: Find a Gilbert's realization of the dynamical structure function. 

First, using the results in [4, Lemma 1], we find a minimal realization (A\, J3 1? Ci, D\) of 
[sQ sP]. When [sQ sP] has I simple poles, using Gilbert's realization fl3| gives 

1 K 

[sQ sP] = V V + lim l s Q sP l 

L ' S — An S^OO 

1=1 % 

where K { = lim s ^ A .(5 — Xi)[sQ sP] and has rank 1 since we are assuming that [sQ sP] has 
simple poles. 

Consider a matrix decomposition of Ki of the following form: 

Ki = EiFu Vi, 



where E % e R p and F % = (Ej E^Ej K x . Then A x = diag{A,} e R lxl , B x 



F'{ F? ... FT 



Ei E 2 ... Ei 



and Di lim^oo [sQ sP]. 
Step 2: Find the maximal number of cancelled poles. 

We define $ as a largest subset of {£?!,-•• , E{\ such that all the elements in $ are mutually 
orthogonal. We also define cj) as the cardinality of $. Computationally, cj) can be obtained using 
the algorithm presented in the Appendix. We claim that cp is equal to the maximum number of 
poles we can eliminate (the proof is in the Appendix). Therefore, the minimal order of [W V] 
is 



2 iV with corresponding minimal realization (A2, B2, C2, /) is restricted to the set of matrices of the form (si — R*)s 1 with 
a constant R* from Proposition |5j 
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As a consequence, the order of the minimal reconstruction is the dimension of A n (the constant 
p) plus the minimal dimension of A 2 2 (obtained above): p + I — 0. 
Step 3: Construct i?* to obtain the minimal reconstruction. 

Once we have <£>, using eq. ( [T3] ) and R* = si — sN, we know that N(\i)\j,j] = implies 
R*\j,j] = A^. Consequently, each element in the set $ will determine at least one element in 
R*. This last fact can be used to construct R* element by element. Once R* is found, we can 
obtain A and B using eq. ( fTO] ). 



IV. Illustrative example 
Example 2: Consider a dynamical structure function [Q,P]: 
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We first compute the McMillan degree of the corresponding transfer function: deg{G} = deg{(7— 
Q)~ l P)} = 4, meaning that a 4 th order state-space model is enough to realize the transfer 
function. It is interesting to look at the minimal order realization consistent with the dynam- 
ical structure function. The different steps of the algorithm proposed in the previous section 
successively yield the following: 

Step 1: A minimal Gilbert realization of s[Q, P] is 

A 1 = diag{-l, -2, -3, -4}, B 1 = diag{2, 2, 2, 4}, 
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Step 2: By definition, Ei 
Thus, 



C\V{ where Vi e M 4 has 1 in its i th position and zero otherwise 
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Furthermore, is 1 and the order of a minimal realization of the given dynamical structure 
function is p + l — = 3 + 4 — 1 = 6. Hence, the system must contain at least 3 hidden states. 
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Step 3: R* can be chosen as diagja, —1, —1}, diag{— 2, a, —2}, diag{— 3, —3, a}, or diag{— 4, —4 
for any a G M. 

The reconstructed networks are represented in Fig. [3] There are three measured (red) nodes, 
labeled 1,2,3 and by the analysis above, there are at least three hidden nodes such that the 
corresponding realization is consistent with the dynamical structure function. The red connections 
between measured nodes are the same for all candidate networks which is in accordance with 
Proposition [TJ Dashed lines correspond to the connections between hidden and measured nodes. 




Fig. 3. Topologies corresponding to the four [Q, P] minimal realizations. The measured nodes are colored red, while the hidden 
ones blue. Red connections between measured nodes are the same for all the networks due to Proposition [T] Each node has a 
self-loop but we omit it for simplicity. 

From a biological perspective, this indicates that there are at least 3 unmeasured species 
interacting with the measured species. Of course, the "true" biological system might be even 
more complicated, i.e., it might have more than 6 species. Yet, when more states are measured, the 
dynamical structure functions can be easily updated and a new search for a minimal realization 
of the updated system can be performed to reveal the corresponding minimal number of hidden 
states. 

V. Conclusions 

In this paper, we have presented a method for obtaining a minimal order realization consistent 
with a given dynamical structure function. We show that the minimal order realization of a given 
dynamical structure function can be achieved by choosing a constant diagonal matrix R*. This 



September 19, 2012 



DRAFT 



14 



provides a way to estimate the complexity of the system by determining the minimal number 
of hidden states that needs to be considered in the reconstructed network. For example, in the 
context of reconstruction of biological networks from data, it helps to understand the minimal 
number of unmeasured molecules in a particular pathway. 
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Appendix 



Proof of the claim in Step 2 of the proposed algorithm: 

Proof: Using results from Section 4 of fl2| ], if a pole of [sQ sP], say A*, is cancelled by 
N = (sI-R^)s- 1 = C 2 (A 2 -sI)- 1 B 2 +I, then the realization of the cascade (sI-fys'^sQ sP] 
loses observability. In this case, it follows that there exists a nonzero vector wi = [w^^w^^ 
such that 

A x - XJ 
B 2 C X A 2 - XJ 

C\ c 2 

The first equation shows that wij is an eigenvector of A\ corresponding to A^. Since A\ is 



0. 



diagonal, w[. 







L 











six/ 



Therefore, we have 



A 2 -\ i I B 2 
C 2 I 
Noticing that C\W\^ = Ei and that 



W 2 ,i 



= 0. 



A 2 -sI B 2 
C 2 I 



I 

-C 2 (A 2 -sI)- 1 I 

we obtain, since \ ^ is not a pole of N, 

N(\i)Ei = 0. 



A 2 - si B 2 
N{s) 



(13) 
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In summary, designing R* to cancel any pole of [sQ sP] is equivalent to imposing that 
eq. ( fT3] ) holds. The next question is: given [sQ sP] what is the maximal number of poles that 



can be cancelled by N, i.e., what is the largest number of poles for which eq. fll3|) is satisfied? 

To answer this, notice that Ei[j] being nonzero for some j 9 implies that there exists at least 
one nonzero element in the j th row of E{. In this case, satisfying eq. ([13]) imposes that the j th 
diagonal element of N(Xi) is 0, i.e., the j th diagonal element of R* is A*. In other words, a 
nonzero element in Ei corresponds to a fixed value in the corresponding diagonal position in 
R*. Since R* is a constant diagonal matrix then any pair of orthogonal vectors in {E\, • • • , E{\ 
does not intervene in the choice of an element on the diagonal of if. ■ 
Algorithm to find (j) and $: 

As is presented in p6| , an undirected graph is denoted by Q = (V, £) where V = {^i, 
is the set of nodes and £ C V x V is the set of edges. 

For our purposes, we construct an undirected graph Q a using the following rules: 

• A node is associated with each vector in the set {E^ • • • , £7}. There are thus I nodes in 
the considered graph. 

• An undirected edge is drawn between node i and node j if the equality EjEj = is 
satisfied. 

It is easy to see that the maximum cardinality of the set $ corresponds to the maximum 
number of nodes in a complete subgraph K n of the graph Q a . 

Although the problem of finding a largest complete subgraph in an undirected graph is a 
NP-hard problem, methods to this end have been well- studied in [17] ^To our best knowledge, 
for an arbitrary graph, the fastest algorithm has a complexity of 0{2 n ^) [18]. Therefore, we 
can use these methods to obtain a largest complete subgraph and consequently compute the 
corresponding set $ and its corresponding cardinality (j). 



3 Some corresponding MATLAB code can be downloaded from http://www.mathworks.com/madabcentral/fileexchange/19889. 
September 19, 2012 DRAFT 



